% March, 2012
%
% Computes the scattering of a plane wave with eight regularly
% spaced objects
% kh    Helmholtz parameter
% angle Angle of the incoming wave

function eightobj(kh,theta_in)

check_sol = true;
verif2d   = false;
tikzoutput = false;

% number of points per wavelength
ppl = 20;
acc = 1e-10;

nobj            = 4;
% The number of points on each contour.
n               = 100;  
ntot            = n*nobj;

% Grid size for 2d visualization
ngrid           = 50;

contours        = cell(nobj,1);
Ctot            = [];
hmax_phy        = 0;
centers  = [[0,0];
                [5,0];
                [0,5];
                [5,5];
               ];
nobjsx = 2;
nobjsy = 2;

figure(1); hold all
for i=1:nobj
  angle       = 2*pi*rand();
  k           = floor(6*rand());
  [C,len,xin,tmp] = get_geometry_star(n,centers(i,:).',k,angle);
  h_phy       = max(sqrt(C(2,:).^2 + C(5,:).^2));    
  hmax_phy    = max(hmax_phy,h_phy);
  contours{i} = C;
  xc(i)       = mean(C(1,:));
  yc(i)       = mean(C(4,:));
  R(i)        = max(sqrt((C(1,:) - xc(i)).^2 + (C(4,:) - yc(i)).^2 ));
  plot(C(1,:) ,C(4,:) ,'k.')
  Ctot        = [Ctot,C];          
end    

xctot       	= mean(Ctot(1,:));
yctot       	= mean(Ctot(4,:));
Rtot            = max(sqrt((Ctot(1,:) - xctot).^2 + (Ctot(4,:) - yctot).^2 ));

lambda = 2*pi/kh;

if verif2d
  figure(2), hold all
  for i=1:nobj
    plot(contours{i}(1,:),contours{i}(4,:),'k')
  end
  figure(3), hold all
  for i=1:nobj
    plot(contours{i}(1,:),contours{i}(4,:),'k')
  end
end

% Check the size of the discretization vs the wavelength
fprintf(1,'h / lambda = %g\n', kh*hmax_phy/(2*pi))

if kh*hmax_phy/(2*pi) > 1./10
  error('mesh too coarse')
end

if verif2d
  xx      = linspace(xctot - 3*Rtot,xctot + 3*Rtot,ngrid);
  yy      = linspace(yctot - 3*Rtot,yctot + 3*Rtot,ngrid);
  [xxx,yyy] = meshgrid(xx,yy);
  %%% Determine which nodes in the visualization mesh are inside the contour
  %%% by evaluating Cauchy integrals.
  %%% Note that one integral is evaluated for each node in the mesh so this
  %%% can be slightly expensive. 
  %%% If you want a finely resolved visualization mesh, it might be worth
  %%% it to precompute "ind" and store it.
  zzz = [reshape(xxx,1,numel(xxx));...
	 reshape(yyy,1,numel(yyy))];
  nz  = size(zzz,2);
  zz =     (ones(nz,1)*Ctot(1,:) - zzz(1,:)' * ones(1,ntot)) + ...
	   1i * (ones(nz,1)*Ctot(4,:) - zzz(2,:)' * ones(1,ntot));
  h   = 2*pi / n;
  dz  = h*Ctot(2,:)' + h*1i*(Ctot(5,:)');
  uu  = (1/(2*pi*1i))*(1./zz)*dz;
  ind = find(abs(uu) > 0.01);
else
  % Verify the solution on a circle
  theta_zzz = linspace(0,2*pi,500).';
  xxx       = xctot + 3*Rtot*cos(theta_zzz);
  yyy       = yctot + 3*Rtot*sin(theta_zzz);
  zzz = [reshape(xxx,1,numel(xxx));...
	 reshape(yyy,1,numel(yyy))];

  figure(1)
  plot(zzz(1,:),zzz(2,:),'k')		
end


% Assemble the matrix
matsA    = cell(nobj,1);
matA   = zeros(ntot);
for i=1:nobj
  i0 = (i-1)*n + 1; i1 = i0 + n - 1;
  matAii = get_A_single_diag(contours{i},1:n,kh);
  matsA{i}          = matAii;
  if (check_sol)
    for j=1:nobj
      i0 = (i-1)*n + 1; i1 = i0 + n - 1;
      j0 = (j-1)*n + 1; j1 = j0 + n - 1;
      if (i == j)
	matA(i0:i1,j0:j1) = matAii; 
      else
	matA(i0:i1,j0:j1) = get_A_offd_noquad(Ctot,i0:i1,j0:j1,kh);
      end
    end
  end
end
level = 1;

% Hierarchical decomposition
level = 1;
res = [[1;ntot]];
anim.isave = 1;
anim.fname = '8obj_step%05d.png';
anim.axis  = [xctot - (Rtot + 1),xctot + (Rtot + 1),yctot - (Rtot + 1),yctot + (Rtot + 1)];
Dplts=cell(2,1);
Csss=cell(2,1);
Contours={contours};
while nobj>=2
  fprintf(1,'\n\n Level %d \n',level)
  level = level + 1;
  acounter=1;
  figure(level); hold all; axis off
  for i=1:nobj
    plot(contours{i}(1,:),contours{i}(4,:),'.')
	if (tikzoutput)
		fid = fopen(['obj4' num2str(level) num2str(i) '.txt'],'wt');
		CCC=[contours{i}(1,:);contours{i}(4,:)];
		fprintf(fid,'%12.8f  %12.8f\n',CCC);
		fclose(fid);
		fid = fopen(['obj4' num2str(level) num2str(i) '.tex'],'wt');
		fprintf(fid,'%d\n',size(contours{i},2));
		fclose(fid);
	end
  end
  [contours_up,matsA_up,nobjsx,nobjsy,Dplt,Css] = factor_level(nobjsx,nobjsy,contours, ...
						  matsA,kh,acc,ppl,level);    
						
  Dplts{level}=Dplt;
  Csss{level}=Css;
  acounter=acounter+1;
  Contours{level}=contours_up;
  contours = contours_up;
  matsA    = matsA_up;
  nobj     = nobjsx*nobjsy;
  
  npts = 0;
  for i=1:nobj
    Ctmp = contours{i};
    plot(Ctmp(1,:) ,Ctmp(4,:) ,'.')
    npts = npts + size(Ctmp,2);
  end
  res = [res,[level;npts]];
  
end

level = level + 1;
figure(level); hold all; axis off
for i=1:nobj
  plot(contours{i}(1,:),contours{i}(4,:),'.')
 % plot(Dplt(1,:),Dplt(2,:),'g.');

  if (tikzoutput)
  	fid = fopen(['obj4' num2str(level) num2str(i) '.txt'],'wt');
  	CCC=[contours{i}(1,:);contours{i}(4,:)];
  	fprintf(fid,'%12.8f  %12.8f\n',CCC);
  	fclose(fid);
  	fid = fopen(['obj4' num2str(level) num2str(i) '.tex'],'wt');
  	fprintf(fid,'%d\n',size(contours{i},2));
  	fclose(fid);
  end
end

%print(gcf, '-dpng', sprintf('8obj_level%02d.png',level-1));

[D,h_D]          = get_circle(contours,1,kh,ppl,1.);
[Is,Cs,V,matP,D] = get_skeleton(kh,contours{1},matsA{1},acc,D,h_D);

level = level + 1;
figure(level); hold all; axis off
for i=1:nobj
  plot(Cs(1,:),Cs(4,:),'.')
  plot(D(1,:),D(2,:),'g.');

  if (tikzoutput)
  	fid = fopen(['obj4' num2str(level) num2str(i) '.txt'],'wt');
  	CCC=[Cs(1,:);Cs(4,:)];
  	fprintf(fid,'%12.8f  %12.8f\n',CCC);
  	fclose(fid);
  	fid = fopen(['obj4' num2str(level) num2str(i) '.tex'],'wt');
  	fprintf(fid,'%d\n',size(Cs,2));
  	fclose(fid);
  end
end
%print(gcf, '-dpng', sprintf('8obj_level%02d.png',level-1));


% Direct solution
v     = -incoming_wave(kh,theta_in,Ctot(1,:),Ctot(4,:)); v = v.';
sigma = matA\v;
phi_scattered_ref  = evalpot(zzz,Ctot,sigma,kh);

% Approximate solution
v     = -incoming_wave(kh,theta_in,Cs(1,:),Cs(4,:)); v = v.';
sigma = matP*v;
phi_scattered      = evalpot(zzz,Cs,sigma,kh);

err = phi_scattered - phi_scattered_ref;

if verif2d
  err      = reshape(err,ngrid,ngrid);
  step     = 4;
  [Inew,idx] = sort(Is);
  Cplot      = Cs(:,idx);
  Cplot      = Cplot(:,1:step:end);
  
  figure(4), hold all;
  err(ind) = NaN;
  contourf(xxx,yyy,log10(abs(err)),-16:2:0)
  cblabel = colorbar('vert');
  plot(Cplot(1,:),Cplot(4,:),'b.','markersize',20)
  plot(D(1,:),D(2,:),'r','linew',2)
  axis equal
  alabel = get(gcf,'CurrentAxes');
  % tlabel = title('outgoing field');
  set(alabel,'FontName','cmr10','FontSize',15,'xscale','lin','yscale','lin','zscale','lin');
  set(cblabel, 'FontName', 'cmr10', 'FontSize', 15);
  % set(tlabel,'FontName', 'cmr10','FontSize',15)
  print(gcf, '-dpng', '2obj_err2.png');

else
  fprintf(1,'\n\n ERROR :  %g \n',max(abs(err)))
end

mycol={'r.','g.','b.','k.'};
% original %%%%%%%%%%%%%%%%%
figure,hold on
axis equal
for ii=1:4
C1=Contours{1}{ii};
plot(C1(1,:),C1(4,:),mycol{ii});
Cout=[C1(1,:);C1(4,:)];
fid = fopen(['obj4o' num2str(ii) '.txt'],'wt');
fprintf(fid,'%12.8f  %12.8f\n',Cout);
fclose(fid);
fid = fopen(['obj4o' num2str(ii) '.tex'],'wt');
fprintf(fid,'%d\n',size(Cout,2));
fclose(fid);
end

% nr 1 %%%%%%%%%%%%%%%%%
figure, hold on
axis equal
%1
C2=Csss{2}{1,1};
plot(C2(1,:),C2(4,:),mycol{1});
Cout=[C2(1,:);C2(4,:)];
fid = fopen(['obj44' num2str(1) '.txt'],'wt');
fprintf(fid,'%12.8f  %12.8f\n',Cout);
fclose(fid);
fid = fopen(['obj44' num2str(1) '.tex'],'wt');
fprintf(fid,'%d\n',size(Cout,2));
fclose(fid);
%2
C2=Csss{2}{2,1};
plot(C2(1,:),C2(4,:),mycol{2);
Cout=[C2(1,:);C2(4,:)];
fid = fopen(['obj44' num2str(2) '.txt'],'wt');
fprintf(fid,'%12.8f  %12.8f\n',Cout);
fclose(fid);
fid = fopen(['obj44' num2str(2) '.tex'],'wt');
fprintf(fid,'%d\n',size(Cout,2));
fclose(fid); ...
%3
C2=Csss{2}{1,2};
plot(C2(1,:),C2(4,:),mycol{3});
Cout=[C2(1,:);C2(4,:)];
fid = fopen(['obj44' num2str(3) '.txt'],'wt');
fprintf(fid,'%12.8f  %12.8f\n',Cout);
fclose(fid);
fid = fopen(['obj44' num2str(3) '.tex'],'wt');
fprintf(fid,'%d\n',size(Cout,2));
fclose(fid); ...
%4
C2=Csss{2}{2,2};
plot(C2(1,:),C2(4,:),mycol{4});
Cout=[C2(1,:);C2(4,:)];
fid = fopen(['obj44' num2str(4) '.txt'],'wt');
fprintf(fid,'%12.8f  %12.8f\n',Cout);
fclose(fid);
fid = fopen(['obj44' num2str(4) '.tex'],'wt');
fprintf(fid,'%d\n',size(Cout,2));
fclose(fid);

D2a=Dplts{2}{1,1};
D2b=Dplts{2}{2,1};
plot(D2a(1,:),D2a(2,:),mycol{1})
plot(D2b(1,:),D2b(2,:),mycol{2})
fid = fopen(['obj44D1.txt'],'wt');
fprintf(fid,'%12.8f  %12.8f\n',D2a);
fclose(fid);
fid = fopen(['obj44D2.txt'],'wt');
fprintf(fid,'%12.8f  %12.8f\n',D2b);
fclose(fid);
D2a=Dplts{2}{1,2};
D2b=Dplts{2}{2,2};
plot(D2a(1,:),D2a(2,:),mycol{3})
plot(D2b(1,:),D2b(2,:),mycol{4})
fid = fopen(['obj44D3.txt'],'wt');
fprintf(fid,'%12.8f  %12.8f\n',D2a);
fclose(fid);
fid = fopen(['obj44D4.txt'],'wt');
fprintf(fid,'%12.8f  %12.8f\n',D2b);
fclose(fid);

% nr 2 %%%%%%%%%%%%%%%%%
figure, hold on
axis equal
for ii=1:2
C3=Csss{3}{ii};
plot(C3(1,:),C3(4,:),mycol{ii});
Cout=[C3(1,:);C3(4,:)];
fid = fopen(['obj42' num2str(ii) '.txt'],'wt');
fprintf(fid,'%12.8f  %12.8f\n',Cout);
fclose(fid);
fid = fopen(['obj42' num2str(ii) '.tex'],'wt');
fprintf(fid,'%d\n',size(Cout,2));
fclose(fid);
end
for jj=1:2
	D3=Dplts{3}{jj};
	plot(D3(1,:),D3(2,:),mycol{jj})
	fid = fopen(['obj42D' num2str(jj) '.txt'],'wt');
	fprintf(fid,'%12.8f  %12.8f\n',D3);
	fclose(fid);
end

% nr 3 %%%%%%%%%%%%%%%%%
figure,hold on
axis equal
plot(Cs(1,:),Cs(4,:),mycol{1})
plot(D(1,:),D(2,:),mycol{1})
Cout=[Cs(1,:);Cs(4,:)];
fid = fopen(['obj411.txt'],'wt');
fprintf(fid,'%12.8f  %12.8f\n',Cout);
fclose(fid);
fid = fopen(['obj411.tex'],'wt');
fprintf(fid,'%d\n',size(Cout,2));
fclose(fid);
fid = fopen(['obj41D1.txt'],'wt');
fprintf(fid,'%12.8f  %12.8f\n',D);
fclose(fid);
keyboard

return
